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ABSTRACT 

Aims. We clarify the response of extrasolar planetary systems in a 2: 1 mean motion commensurability with masses ranging from the 
super Jovian range to the terrestrial range to stochastic forcing that could result from protoplanetary disk turbulence. The behaviour 
of the different libration modes for a wide range of system parameters and stochastic forcing magnitudes is investigated. The growth 
of libration amplitudes is parameterized as a function of the relevant physical parameters. The results are applied to provide an expla- 
nation of the configuration of the HD 1283 11 system. 

Methods. We first develop an analytic model from first principles without making the assumption that both eccentricities are small. 
We also perform numerical N-body simulations with additional stochastic forcing terms to represent the effects of putative disk tur- 
bulence. 

Results. We isolate two distinct libration modes for the resonant angles. These react to stochastic forcing in a different way and be- 
come coupled when the libration amplitudes are large. Systems are quickly destabilized by large magnitudes of stochastic forcing but 
some stability is imparted should systems undergo a net orbital migration. The slow mode, which mostly corresponds to motion of the 
angle between the apsidal lines of the two planets, is converted to circulation more readily than the fast mode which is associated with 
oscillations of the semi-major axes. This mode is also vulnerable to the attainment of small eccentricities which causes oscillations 
between periods of libration and circulation. 

Conclusions. Stochastic forcing due to disk turbulence may have played a role in shaping the configurations of observed systems in 
mean motion resonance. It naturally provides a mechanism for accounting for the HD 12831 1 system for which the fast mode librates 
and the slow mode is apparently near the borderline between libration and circulation. 

Key words. Turbulence - Celestial mechanics - Planetary systems: formation - Planetary systems: protoplanetary disks - GJ876: 
planetary systems - HD128311: planetary systems: formation 

1 . Introduction sue (eg. i Ganrnild [19961) . Because the level of MRI turbulence 

and associated density fluctuations are very uncertain, we con- 

Of the recendy discovered 335 extrasolar planets, at least 75 are ^^^^^ stochastic force amplitudes ranging over several orders of 

in multiple planet systems (Schneider 2009). About 10% of these magnitude 
are in or very close to a resonant configuration where two planets 

show a mean motion commensur abilit y, with four sys tems in or The influence of stochastic forces resulting from the grav- 

near a 2:1 resonance dUdrv et al.ll2007HReipurfli et al.ll2007l) . itational field produced by the density fluctuations associated 

Resonant configurations can be established by dissipative with MRI turbulence on migrating planets was explored first by 

forces acting on the planets whic h lea d to convergent migra - iNelson & Papaloiz ou (2004). They considered MRI simulations 

tion (see for exampIe lLee & Fe^Mi^^^dto^^tFinm "^''^f'^y "^f the result that the simulation ran only for a re a- 

An anomalous effective kinematic viscosity ly of order ~ IQ-^, lively small number of orbits. To consider much longer evolu- 

but with considerable uncertainty, has been inferred from obser- tion.times, as done here, a simpler parameterized model of the 

vations of accretion rates onto the central protostars. The mag- orcmg is nee e . 

neto rotational instability (MRI) is thought to be responsible for Re cent studies by lAdams et alj (l2008l) and iLecoanet et al.l 

this anomalous value of v which is often characterized using the ( l2008h have been made in order to estimate the lifetimes of the 

Shakura & Sunyaev (1973) a — v/{csVl) parameter, with Cs be- mean motion resonances in the GJ876 system when it is per- 

ing the local sound speed and 2tt/VI being the orbital period. For turbed by a sequence of stochastic kicks. They found that res- 

the most likely situation of small or zero net magnetic flux, re- onances were disrupted within expected disk lifetimes for suf- 

cent MHD simulations have indicated a w 10^^ — 10^"^. Even ficiently large forcing. They assume a fixed orbit for the outer 

assuming adequate ionization, this value is uncertain when real- planet and gave an interesting discussion of the interactions of 

istic values of the actual transport coefficients are employed be- the two planets based on a pendulum equation with an additional 

cause of numerical resoluti on issues (see Fromang & Papaloizou ad hoc stochastic forcing term. The aim of the work presented in 

120071: iFromang et alj2007h . Which parts of protoplanetary disks this paper is to provide a more complete study of the full set of 

are adequately ionized, or consitute a dead zone, is also an is- celestial mechanics equations and to describe the interplay be- 
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tween apsidal corotation and mean motion resonance as well as 
to consider a wide range of planet masses and stochastic force 
amplitudes. In addition, we derive a prescription for incorporat- 
ing continuous stochastic forcing terms that allows for a general 
autocorrelation function with associated correlation time. Their 
magnitudes are related to properties of the protoplanetary disk 
and physical scaling laws are found. This is of particular im- 
portance as the factors that control the strength of the stochas- 
tic forcing are not well constrained. We also use our formalism 
to estimate the stochastic diffusion rates of the orbital elements 
from first principles. 

We further remark that Adams et al. (2008) gave a discus- 
sion of the diffusion of reson ant angles that do not satisfy the 
d'Alambert condition (see eg. lHamiltonlll994h which is equiva- 
lent to the requirement of rotational invariance. Thus their ap- 
pearrance cannot be straightforwardly connected to the basic 
equations or angles discussed in this paper unless one assumes 
that the variation of the longtitude of pericentre of one or both 
of the interacting planets can be neglected. However, the system 
lifetimes derived from their numerical work is broadly consistent 
with ours for the parameter regime they considered. 

We find that stochastic forcing readily produces systems in 
mean motion resonance with broken apsidal corotation. An ad- 
ditional aim of this paper is to use this feature to construct sce- 
narios involving convergent migration and stochastic forcing to 
account for the HD12831 1 system. 

The plan of the paper is as follows. In section 2 we present 
the basic equations governing two interacting planets subject to 
external stochastic forces. We specialize to the case where the 
planets are either in or near a mean motion commensurability 
and so retain only the two angle variables that vary on a time 
scale much longer than the orbital one. Considering the case 
where these angles undergo small amplitude librations, we iden- 
tify fast and slow libration modes, the former being associated 
with variations in the semi-major axes and the latter with the an- 
gle between the apsidal lines of the two planets. We go on to 
consider the effects of stochastic forcing arising from some ex- 
ternal process such as disk turbulence in section 2.3 deriving the 
diffusion rates for the orbital elements of a single planet and the 
growth rate of the libration amplitudes in the two planet case. 

In section 3 we discuss the origin and numerical implemen- 
tation of the stochastic forcing and their operation in the sin- 
gle planet case. In section 4 we go on to consider the stochastic 
forcing in the two planet case. We consider the conversion of 
the fast and slow modes from libration to circulation for model 
systems of varying total mass ratio and initial eccentricities in- 
cluding GJ876. The attainment of small eccentricities and cou- 
pling to the fast mode results in conversion of the slow mode to 
circulation before the fast mode. Investigations are carried out 
for stochastic diffusion rates, proportional to the mean square 
stochastic force amplitude ranging over several orders of magni- 
tude. The life time of the resonant angles is found to be inversely 
proportional to the diffusion rate except in the case of systems 
with low total mass in the earth mass range. 

In section 5 we exploit the tendency of the slow mode, re- 
lated to the angle between the apsidal lines, to be driven to cir- 
culate while the fast mode still librates in stochastically forced 
systems, to make a model for the formation of the HD128311 
system which may be in such a state and was not readily un- 
derstood in terms of convergent migration models for producing 
the commensurability. We combine the effects of such migra- 
tion and stochastic forcing, showing that during the migration 
phase, while the librations tend to be stabilized, the slow mode 
is readily converted to circulation while the fast mode contin- 



ues to librate. Good agreement is obtained with the somewhat 
uncertain observed orbital configuration. Finally in section 6 we 
summarize and discuss our results. 



2. Basic Equations 

We begin by writing down the equations of motion for a single 
planet moving in a fixed pl ane under a general Hamiltonian H 
in the form (see e.g. ,Snellgrove et aLll2001; .Papaloizou.20 03) 



E 
G 
X 



dH 

dH 
9A 



dH 
dm 



dH 

'dL 
dH 

'dL' 



dH 
'~dE 



(1) 
(2) 

(3) 
(4) 



Here the angular momentum of the planet is G and the energy is 
E. For orbital motion around a central point mass A/ we have 



G = myJgMa{l - e^) and 



E = -'- 



2a 



(5) 
(6) 



where Q is the gravitational constant, a the semi-major axis and 
e the eccentricity. The mean longtitude is A = n{t — to) + nj, 
where n the mean motion, with to being the time of periastron 
passage and vu being the longitude of periastron. 

2. 1. Additional Forcing of a Single Planet 

In order to study the phenomena such as stochastic forcing we 
need to consider the effects of an additional external force per 
unit mass F which may not be described using a Hamiltonian 
formalism. However, as may be seen by considering general co- 
ordinate transformations starting from a Cartesian representa- 
tion, the equations of motion are linear in the components of F. 
Because of this we may determine them by considering forces of 
the form F = {F^, Fy) for which the Cartesian components are 
constant. Having done this we may then suppose that these vary 
with coordinates and time in an arbitrary manner Following this 
procedure we note that when F, as in the above form is constant, 
we may derive the equations of motion by replacing the original 
Hamiltonian with a new Hamiltonian defined through 



H 



H - m{F^, X + Fy y) = H - m{r -F) . 



(7) 



The additional terms proportional to the components of F 
correspond to the Gaussian form of the equations of motion 
jBrouwer & Clemencelll96ll) . 

The various derivatives involving r can be calculated by el- 
ementary means and expressed in terms of E, G, A and vu. One 
thus finds additional contributions to the equations of motion 
([T]i - (|4|i, indicated with a subscript F, in the form 



Gp ~ m { 
Ef 



d_ 
dX 
d 



TO (r X F) • 



TO (v • F) 



ZUp 



mn— (r • F) 

dA 

d 

— TO^^ (r • F) or equivalently 



(8) 

(9) 
(10) 
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Wp 



F„ 1 



1 



1 — a 



sin f — Fr cos / 



A. = -m(^- + n-j (r.F) 



(11) 
(12) 
(13) 



where the true anomaly / is defined as the difference between 
the true longitude and the longitude of periastron, f = 9 — zu. 
Note that from (|9|l we obtain 



2ah _ 2(F^esin/ + Fe(l + ecos/)) 
3ri nVT 



nv i — 

and from dHJ together with (|9]l we obtain 



G{2EG + GE) 



(14) 



(15) 



In the limit e <C 1 this becomes (ignoring terms 0(e) and 
smaller) 

ep = Fr— sin f + Fg— 2 cos f. (16) 

an an 

Furthermore in this limit we may replace / by / = A — to = 

n{t - to). 

We remark that the above formalism results in equation (fTTT l 
which gives an expression for zu p that indicates that this quantity 
diverges for small e as 1 /e. We comment that, as is well known, 
this aspect results from the choice of coordinates used and is 
not associated with any actual singularity or instability in the 
system. This is readily seen if one uses h = esintu, and k = 
e cos zu as dynamical variables rather than e and zu. The former 
set behave like Cartesian coordinates, while the latter set are the 
corresponding cylindrical polar coordinates. When the former 
set are used, potentially divergent terms cx 1/e do not appear. 
This can be seen from (fTTT) and (fT6b which give in the small e 
limit 



1 



1 



hp = —Fr — sinA + Fe — 2 cos A 
an an 

kp = Fr — sinA + i^e — 2 cos A. 

an an 



(17) 
(18) 



Abrupt changes to zu may occur when h and k pass through the 
origin in the {h, k) plane. But this is clearly just due to a coordi- 
nate singularity rather than a problem with the physical system 
which changes smoothly as this transition occurs. The abrupt 
changes to the zu coordinate occur because very small perturba- 
tions to very nearly circular orbits produce large changes to this 
angle. 

2.2. Multiple Planets 

Up to now we have considered a single planet. However, it is a 
simple matter to generalize the above formalism so that it ap- 
plies to a system of two planet s. He re we follow closel y t he d is- 
cussions in iPapaloizoul (l2003h and [Papaloizou & Szuszkiewi^ 
t2005J . Excluding stochastic forcing for the time being, we start 
from the Hamiltonian formalism des cribing their m utual inter- 
actions using Jacobi coordinates (see lSinclaiij|1975h . In this for- 
malism the radius vector r2 , of the inner planet of reduced mass 
m2 is measured from M and that of the outer planet, ri, of re- 
duced mass nil is referred to the centre of mass of M and 7712. 



Thus from now on we consistently adopt a subscripts 1 and 2 for 
coordinates related to the outer and inner planets respectively. 

The required Hamiltonian correct to second order in the 
planetary masses is given by 



^ Qmim2 



Qmim2Yi 



GMimi gM2m2 



1-2 



TC2 



|ri| 



(19) 



Here Mi = M + mi, M2 = M + 7712 and ri2 = 
r2 — ri. The Hamiltonian can be expressed in terms of 
Ei,Gi,zUi, Xi,i — 1,2 and the time t. The energies are given 
by Ei = —QmiMi/{2ai), and the angular momenta Gi — 
nniyjQMiai{l — e|) with a.; and denoting the semi-major 
axes and eccentricities respectively. The mean motions are ni = 

^/gM~^. 

The Hamiltonian may quite generally be expanded in a Fourier 
series involving linear combinations of t he three angular differ- 
ences Ai — vji, i = 1,2 and zui — zu2 (eg. lBrouwer & Clemencd 
[19611) . 

Near a first order p + 1 : p resonance, we expect that both 
(f>i = {p + l)Ai — pX2 — VJ2, and (j)2 = [p + l)\i ~ p\2 ~ 
zui, will be slowly varying. Following standard prac tice (see eg. 
|Papaloizoull2003l: iPapaloizou & Szuszkiewiczll2005h only terms 
in the Fourier expansion involving linear combinations of (f>i and 
(j>2 as argument are retained because only these are expected to 
lead to large long term perturbations. 

The resulting Hamiltonian may then be written in the general 
form H = El + E2 + H12, where 



^12 — 2^ (^fe,i 



ai 



Oi 
02 



,61,62 cos{l(t)i+k(j)2), (20) 



where in the above and similar summations below, the sum 
ranges over all positive and negative integers (fc, /) and the di- 
mensionless coefficients Gk,i depend on 61,62 and the ratio 
ai/O'2 only. We also replace Mi by M. 

2.2.1 . Equations of Motion 

The equations of motion for each planet can now be 
easily derived that take into account the contributions 
due to their to m utual interactions (see 'Papaloizou' "2003t 
Papa loizou & Szus zkiewicz 2005) and contributions from dHJ - 
( fT3T l. The latter terms arising from external forcing are indicated 
with a subscript F. We thus obtain to lowest order in the perturb- 
ing masses. 



dni 
~dt' 



dn2 



dei 
~dt 



3{p + l)n{m,2 
M 
dni \ 
~dt)p 
?tpn^mia2 



'^Ckjik + I) sin{l(f)i + k(j)2 



Mai 
dn2 
~dt 



"^CkAk + I) smil(f)i + k(t>2) 



(21) 



(22) 



Ckd sm{l<j)i + k(j)2) (23) 



»7i2"-i yl-ef 
eiM 



k-{p+l)ik + l) [l-Jl-ej 



' dei 
IT 
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de2_ 
~dt 



-T-r ■ > , C^fc.i sin(i0i + k(t>2) 



= {p+l)ni~pn2-^{Dkj+Ek,i)cos{l(j)i 



dt 



Here 



dt 



2(p + l)nia^TO2 d 



(24) 

42) 
(25) 

(26) 



(Cfe,,/ai) 



and 



M 502^^'=''/"^)' 
"1^2 ((p + 1)(1 - el) - py/l - QQ 



(27) 



k,l 



pn2a2mi 



eiM 



9ei 



2 dCi 



k,l 



0162*1 9e2 

(p + l)niTO2 (1 - 6? - - ei) dCk.i 
61 M c)ei 
n2a2mi {{p+l)^/l - el - p{l - 6^)) ^(^^^ 



(28) 



(29) 



0162^/ 9e2 

Note that 02 — 0i = 1372 — nji = Aiu is the angle between 
the two apsidal lines of the two planets. We also comment that, 
up to now, we have not assumed that the eccentricities are small 
and that, in additional to stochastic contributions, the external 
forcing terms may in general contain contributions from very 
slowly varying disk tides but we shall not consider these further 
in this article. 

2.2.2. Modes of Libration 

We first consider two planets in resonance with no external 
forces acting in order to identify the possible libration modes. 
We then consider the effects of the addition of external stochas- 
tic forcing. In the absence of external forces equations (1211 1 - ( |26] l 
can have a solution for which and Ci are constants and the an- 
gles (pi and (/)2 are zero. In general other values for the angles 
may be possible but such cases do not occur for the numerical 
examples presented below. When the angles are zero equations 
(I25]) and (l26b pro vide a relationship between 61 and 62 (see eg. 
lPapaloizoul2003h . 

We go on to consider small amplitude oscillations or libra- 
tions of the angles about their above equilibrium state. Because 
two planets are involved there are two modes of oscillation 
which we find convenient to separate and describe as fast and 
slow modes. Assuming the planets have comparable masses, the 
fast mode has libration frequency cx ^/m and the slow mode has 
libration frequency cx m. These modes clearly separate as the 
planet masses are decreased while maintaining fixed eccentrici- 
ties. 



2.2.3. Fast Mode 

To obtain the fast mode we linearize (I2TI 1 - ( l26b and neglect sec- 
ond order terms in the planet masses. This is equivalent to ne- 
glecting the variation of Dkj,Ekj and Fkj in equations ( |25] ) 
and ( l26b which then require that (pi = (/)2 very nearly for this 
mode. Noting that for linear modes of the type considered here 
and in the next section, equations (l2ni-(l24li imply that rii and ii 
are proportional to linear combinations of the librating angles, 
differentiation of either of equations ( |25] l or ( |26] l with respect to 
time then gives for small amplitude oscillations 



dt^ 
where 



'If 



+ uff,p, = 0, (i = 1,2), 



(30) 



3p n2m2 
M 



0.2^1 

aim2 



Y,CkAk + l? (31) 



and we have used the resonance condition that (p+ 1)711 ~ P2i^2 
which is satisfied to within a correction of order a/ mi/M. Note 
that for this mode the fact that (pi = (p2 very nearly, implies that 
•072 — ci7i is small. Thus that quantity does not participate in the 
oscillation. 



2.2.4. Slow Mode 

In this case we look for low frequency librations with frequency 
oc mi. Equations (|25T l and ( |26] l imply that, for such oscillations, 
to within a small relative error of order to 1 /Af , (p + 1 ) n 1 — pn2 
throughout. 

Equations (1211 1 and (l22l i then imply that the two angles are 
related by 02 = l3(pi , 

where /3 = -J^CkA^ + i)l)/{J2CkAk + l)k). Subtracting 
equation ( |26] l from equation dZSl ), differentiating with respect to 
time and using this condition results in an equation 

for Q = (p2 ^ <Pi — '^2 — '^1 — An7 



where 



0, 



^iiTi2y/l — e\ dW n2a2miy/l — e\ dW 
eiM{l-l3) dei'^^^ ai62M(l - /3) 'de^ 



(32) 



(33) 



with 



ai = 



Ct2 = 



and 



W = 



k-{p+l){k + l) [l-Jl-ei 



Ck,. 



l+p{k + l) 



(kp + l) 



I nia2m2y/T' 
y ai62M 



n2miy/l - el d 



de2 



eiM 



dei 



Although the expressions for the mode frequencies are compli- 
cated, the fast frequency scales as the square root of the planet 
mass and the slow frequency scales as the planet mass indepen- 
dent of the magnitude of the eccentricities while both scale as the 
characteristic mean motion of the system. Furthermore although 
we have considered small amplitude librations and accordingly 



H. Rein and J.C.B. Papaloizou: Mean Motion Resonances and Stochastic Forcing 



5 



obtained harmonic oscillator equations, the treatment can be ex- 
tended to consider finite amplitude oscilations and generalized 
pendulum equations as long as the two mode frequencies can 
be separated. However, we shall not consider this aspect further 
here. 



2.2.5. Librations with External Forcing 

When external forcing is included source terms appear on the 
right hand sides of equations ( |30] | and (l32T i. We shall assume 
that the forcing terms are small so that terms involving products 
of these and both the libration amplitudes and the planet masses 
may be neglected. Then in the case of the slow mode, repeating 
the derivation given above including the forcing terms, we find 
that (|32|i becomes 



dt 



{W2F — Ti'ip) 



(34) 



The quantities ujiF are readily obtained for each planet from 
( fTTT i. From this we see that for small eccentricities, vuip oc l/e^, 
indicating large effects when is small. As already discussed 
in section 12.11 this feature arises from a coordinate singularity 
rather than physically significant changes to the system. 

A similar description may be found for the fast mode. In this 
case, neglecting terms of the order of the square of the planet 
masses or higher, one may use equations (|25l l and (|26] l to obtain 
an equation for Q = 0i in the form 



d^Q 
dt^ 



'If- 



d 
di 

+ {p+ l)niF -pn2F- 



(35) 



If equations for quantities that are regular functions of 
the Cartesian like coordinates {hi = sin t<7i and fc^ = 
eicoswi, i — 1,2) are found, the source terms arising from 
the external forcing are regular as 0. To illustrate this we 

consider the pair 

y = eie2 sin = ft,2fci — /iifc2 and 

Z = 6162 cos ^ = /Ci/C2 ~ /il^2- 

It is readily seen that their time derivatives are given by 

. _ ^(6162) 



dt 



cos( 



■ sm( + 6162 cosC 



dc 

dt 



(36) 



d(eie2) 



dt 



e2eiF + eie2F 



d( 

6162 smQ— (37) 
dt 



Here \u denotes evaluation without external forcing. There is 
now no divergence of the external forcing terms as — > 0. If as 
above small amplitude slow mode librations of C, about zero are 
considered, given that without external forcing cx the first 
term on the right hand side of (|36] | is either quadratic in (, or 
proportional to the product of the exernal forcing and libration 
amplitudes and so, as above may be neglected. The same applies 
to last term on the right hand side of ( |37] ). In addition cos C, may 
be replaced by unity. On taking the time derivatives of ( l36b and 
JJTI i, we obtain under the same approximation scheme 



ei62 



dt^ 



— [6162 {VJ2F - ^1f)] 



62 



d^ei 



dt^ 



ei 



^^62 



dt^ 



— [eie2F + 6261^] 



(38) 
(39) 



Making use of ( |32] | and its counterparts for e'l and 62 for the 
unforced motion we obtain 



y 



dz 



''5z 



d 

'dt 
d_ 

^ It 



[6162 (ti72F - t^if)] 



\e1e2F + e2eiF 



(40) 



(41) 



where 6z = z — zq, zq being the value of 6162 at the centre of 
the oscillation in z. These may be used as an alternative to ( l34l i. 
and are without potentially singular forcing terms. We note that 
for the purpose of this paper, both descriptions provide exactly 
the same physical content. 

Equations ( |35] | and ( |34] | (or ( |40] | together with dTO ) form 
a pair of equations for the stochastically forced fast and slow 
modes respectively. We comment that this mode separation is 
not precise. However, it can be made so by choosing appropri- 
ate linear combinations of the above modes. Numerical results 
confirm that Q predominantly manifests the fast mode and C or 
y the slow mode, so we do not expect such a change of basis to 
significantly affect conclusions. 

We further comment that because (piF contains d72F but not 
tijiF, for small eccentricities there are only potential forcing 
terms oc l/e2 that occur when forcing is applied to the inner 
planet. As this planet has the larger ecentricity for the situations 
we consider, small eccentricities are not found to play any sig- 
nificant role in this case. 

Each mode responds as a forced oscillator We suppose the 
forcing contains a stochastic component which tends to excite 
the respective oscillator and ultimately convert libration into cir- 
culation. But we stress that the above formulation as well as 
developments below assume small librations, so we may only 
assess the initial growth of oscillation amplitude. However, in- 
ferences based on the structure of the non linear governing equa- 
tions and an extrapolation of the linear results enable successful 
comparison with numerical results. 

2.3. Stochastic Forces 

We assume that turbulence causes the external force per unit 
mass {Fr, Fg) acting on each planet to be stochastic. For sim- 
plicity we shall adopt the simplest possible model. Regarding 
the components of the force, per unit mass expressed in cylindri- 
cal coordinates to be a function of time t, we assume any one of 
these satisfies the relation Fi{t)Fi{t') = l^Ff) g{\t - t'\) where 
the autocorrelation function g{x) is such that g{x)dx — Tc, 

where Tc is the coiTelation time and \/{Ff) is the root mean 
square value of the i component. It should be noted that an en- 
semble average is implied on the left hand side. Again for sim- 
plicity we assume that different components acting on the same 
planet as well as components acting on different planets are un- 
coiTelated. We note that in general the root mean square values 
as well as may depend on t, but we shall not take this into ac- 
count here and simply assume that these quantities are constant 
and the same for each force component. 

We note the stochastic forces make quantities they act on 
undergo a random walk. Thus if for example A — Fi for some 
quantity A (note that constants or slowly varying quantities orig- 
inally multiplying Fi may be absorbed by a redefinition of A 
and so do not materially affect the discussion given below), the 
square of the change of A occuring after a time interval t is given 
by 



F,{t')F,{t")dt'dt" 
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"'0 



{Fl)g{\t' ^t"\)dt'dt" ^ D,t, 



(42) 



Here we take the limit where t/Tcis very large corresponding to 
an integration time of very many correlation times and 
Di = 2 (-Ff ) Tc is the diffusion coefficient. 

When the evolution of a stochastically forced planetary orbit 
is considered, it is more appropriate to consider a model govern- 
ing equation for A of the generic form 



A = Fi sm{nt), 



(43) 



where we recall that 27r/n is the orbital period (but note that 
a different value could equally well be considered). We note in 
passing that, by shifting the origin of time, an arbitrary phase 
may be added to the argument of the sin without changing the 
results given below. One readily finds that equation (|42] | is re- 
placed by 



(AA)-' 



JO 

t pt 



"'0 



Fi{t')F,{t") sin(nt') sin{nt")dt'dt" 
{Ff) g{\t' - t"\) sin(?it') sm{nt")dt'dt" 



(44) 



where 7(n) — g{x) cos{nx)dx. 

Note that when tiTc <C 1, corresponding to the correlation time 
being much less than the orbital period, 7^1. For larger Tc, 
7 < 1 gives a reduction factor for the amount of stochastic dif- 
fusion. For example if we adopt an exponential form for the au- 
tocorrelation function such that 



g{\t'-t"\) = cxp 



\t'-t"\ 



we find 
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1 



1 + n'^T^ 



(45) 



and for the purposes of comparison with numerical work we 
shall use this from now on. 



2.3.1 . Stochastic Forcing of an Isolated Planet 

We begin by considering the effect of stochastic forcing on a 
single isolated planet. In this case we may obtain a statistical es- 
timate for the characteristic growth of the orbital parameters as 
a function of time by integrating equations (fTTl i. ( fTSI ), (TT\ 
and ( fTSI ) with respect to time directly. We may then apply the for- 
malism leading to the results expressed in generic form through 
equations (l42l i - ( |45]) to obtain estimates for the stochastic diffu- 
sion of the orbital elements in the limit of small eccentricity in 
the form 



{Aaf 
[Aef 
{Awf 
{Ahf 
{Akf 



2 



= 2.5 



-fDt 



2.5 jDt 



2.5 



2.5 



•yDt 



■jDt 



(46) 
(47) 
(48) 
(49) 
(50) 



We note again that the 1/e^ dependence of (Aro)^ which 
arises from the coordinate singularity discussed in section 12.11 
does not appear for (A/i)^ and (Afc)^. Note that the definition 
of {h, k) imply consistently with the above that 



{Ahf = (Ae)^ sin^ w + {Azof cos^ w 
{Akf = (Ae)^cos^ti7 + e^(Atn)^sin^n7. 



(51) 
(52) 



2.3.2. Stochastic Variation of the Resonant Angles in the 
Two Planet Case 

We now consider the effects of stochastic forcing on the resonant 
angles. The expressions are more complicated than in the previ- 
ous section because there are more variables involved. However, 
we basically follow the formalism outlined in section [23] 

While the libration amplitude is small enough for lineariza- 
tion to be reasonable, the evolution is described by equations 
(|35T l and (O (or (gO)! together with (gB). These may be solved 
by the method of variation of parameters. Assuming the ampli- 
tude is zero at t = 0, the solution of equation ( l40b is given by 

y = sm{ujist) I — ^ -dt 



t^ls 

t 



I SySm{ujist) 
■cos{uJist) ' 



dt. 



(53) 



where Sy = d{eie2{'CJ2F — ■djip))/dt. There are corresponding 
expressions that can be obtained from equations (|4TI) and ( |35] | 
for 6z and Q respectively. 

Equation ( |53] | may be regarded as describing a harmonic os- 
cillator whose amplitude varies in time such that the square of 
the amplitude after a time interval t is given by 



* Sy sm{ujist) 



dt 



Sy COs{uJist) 



dt 



UJls 



The corresponding expression from equation (l4Tl l is 

2 



{ASzf 



Sz s\n{uoist) 



dt 



Sz cos(cj/si) 



(54) 



dt 



(55) 



where Sz ~ d{eie2F + eie2F)/dt. And finally the equation 
obtained from equation i35[ is 



Sq sm{ujift) 

UJlf 



dt 



Sq COs{u!ift) 
LOlf 



dt 



(56) 

where Sq = d{4)iF)ldt + {p + l)niF - pn2F- 

We now evaluate the expectation values of these using the 
formalism of section 12.3! For simplicity and as considered nu- 
merically later we shall specialize to the case when stochas- 
tic forces act only on the outer planet (but see section |6] be- 
low). Taking equation ( |54l ), we perform an integration by parts 
neglecting the end point contribution as these are associated 
with subdominant contributions increasing less rapidly than t for 
large t, to obtain 

{Ayf = (A(eie2)sinC)' 

eie2'coiF sm{ujist)dt 



/ eie2'djip cos{LOist)dt 
Jo 



(57) 
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We then find, working in the limit of a small eccentricity ei (but 
not necessesarily 62), from the corresponding equation for 8z 
that(Afc)2 = (Ay)2. 

In dealing with equation ( l56b we neglect (\)\p in S because 
after integration by parts this leads to a contribution on the or- 
der LOif /n smaller than that derived from hip. Thus we simply 
obtain 



dt 



^if 
t 



hiF COs{uJift) 



dt 



(58) 



/o ^if 

We now follow the procedures outlined in section l23] obtaining 



(Aeiea sinC)^ 

(Ag)^ 

where 



If 

Is 



9Djft 

? 2 



and 



and 



1 



2.3.3. Growth of Libration Amplitudes 



(59) 
(60) 

(61) 
(62) 



Equations (|59] l and ( |60] | express the expected growth of the 
resonant angle libration amplitudes as a function of time. We 
remark that these expressions can be simply related to those 
obtained for a single planet. Thus equations ( |46] | and ( |60] | 
applied to the outer planet imply that 

(AQ)V(Aai)^ = 9ip+irnhf/iAalu;fj). 

As we are interested in the case p — 1, the width of the libration 
zone is ^ aiUif/ni, we see that the time for (AQ)^ to reach 
unity is comparable for the semi-major axis to diffuse through 
the libration zone. 

Similarly, for small amplitude librations about fixed 62 
equation ( |59] l gives almost the same result for ( Aei sin C)^ or 
e2^^(A(5z)^ as that obtained for (Ah)^ or (Afc)^ from equation 
( |49] | applied to the outer planet. This indicates that (, being the 
angle between the apsidal lines of the two planets, diffuses in 
the same way as for an isolated outer planet subject to stochas- 
tic forces. Thus, in the small amplitude regime, the way this 
diffusion occurs would appear to be essentially independent of 
whether the planets are in resonance (but see below). 

An important consequence of equation i5% is the behaviour 
of C, for small ei. The latter quantity was assumed constant in 
the analysis. As implied by the discussion of section IZT] abrupt 
changes to C are expected when {y, z) passes through the origin. 
Then even an initially small amplitude libration is converted to 
circulation. Thus if ei is small then equation ( |59l ) indicates that 
a time t ~ Aa\n\e\/ {bDjsC'^), is required to convert libration 
to circulation. This can be small if ei is small. Even if ei is 
not small initially, it is important to note that it also undergoes 
stochastic diffusion (see equation (|47| i) as well as oscillations 
through its participation in libration. Should ef attain very small 
values through this process, then from ( |59l ) we expect the onset 
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HD128311 b 


1.56 


1.109 


476.8 


0.38 


58" 


c 


3.08 


1.735 


933.1 


0.21 





Table 1. Parameters of the model systems considered. The first 
has orbital elements taken from the three planet fit t o GJ876 with 
orbita l inclination to the plane of the sky, i ~ 50° (iRivera et al.l 
120051). The table entries labelled as LM, SE and E have the same 
parameters as the first entry but the planet masses are scaled 
down by constant factors of 6, 60, and 600. Systems with the 
added label HE have larger orbital eccentricities. The final entry 
is a system with the observed elements of HD1283 1 1 dVogt et al.l 



of a rapid evolution of (. Accordingly the attainment of circula- 
tion for this angle, should be related to the diffusion of ef allow- 
ing very small values of that quantity to be attained, rather than 
the direct excitation of libration amplitude. This is particularly 
the case when ef starts from relatively small values. 

In fact, application of ( |59] l and (|60] | to the numerical exam- 
ples discussed below adopting the initial orbital elements indi- 
cates that the diffusion of ( is significantly smaller than Q unless 
ei starts out with a very small value. This would suggest that Q 
reaches circulation before C,. However, this neglects the coupling 
between the angles that occurs once the libration amplitudes be- 
come significant. It is readily seen that it is not expected that 
Q could circulate while ( remains librating as it was initially. 
One expects to recover standard secular dynamics for ( from the 
governing equations (I2TI 1 - dZSl l. when these are averaged over 
an assumed Q — (pi circulating with constant (j. As a libra- 
tion of the initial form would not occur under those conditions, 
we expect, and find, large excursions or increases in the libration 
amplitude of to be correlated with increases in the libration am- 
plitude of Q. This in turn increases the oscillation amplitude of 
the eccentricity ei, allowing it to approach zero. The consequent 
rapid evolution of C, then enables it to pass to circulation. Thus 
the breaking of resonance is ultimately found to be controlled by 
the excitation of large amplitude librations for Q ~ (pi, which 
induce C, to pass to circulation somewhat before Q itself. 

3. Numerical Simulations 

We have performed numerical simulations of one and two planet 
systems that allow for the incorporation of additional stochastic 
forces with the properties described above. These in turn pro- 
vide a simple prescription for estimating the effects of stochastic 
gravitational forces produced by density fluctuations associated 
with disk turbulence. Results have been obtained using both a 
fifth order Runge Kutta method and the Bulirsch Stoer method 
(Stoer & Bulirsch 2002; Press et al. 1992) with fixed as wefl as 
adaptive timesteps. We have checked that results are converged 
and thus do not actually depend on the integrator used. 
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First we discuss the expected scaling of the stochastic forces 
with the physical parameters of the disk and their implementa- 
tion in the n-body integrations. In order to clarify the physical 
mechanisms involved and to check the analytic predictions for 
stochastic diffusion given by equations ( |46] | - ( |49] l we consider 
simulations of a single planet undergoing stochastic forcing first. 
We then move on to consider two planet commensurable systems 
with and without stochastic forcing. We focus on the way a 2:1 
commensurability, corresponding to p = 1 is disrupted, high- 
lighting the various evolutionary stages a system goes through 
as it evolves from a state with a strong commensurability affect- 
ing the interaction dynamics, to one where the commensurabil- 
ity is completely disrupted and in some cases a strong scattering 
occurs. We consider a range of different planet masses and ec- 
centricities (see table[T]i. 

3.1. Stochastic Forces 

In order to mimic the effects of turbulence, for example pro- 
duced by the MRI, it is necessary to calibrate these forces with 
reference to MHD simulations. As described above, the basic 
parameters characterising the prescription for stochastic forcing 
that we have implemented are the root mean square value of 
the force components per unit mass (in cylindrical coordinates) 
\/{Ff) and the auto correlation time Tc- 

From our analytic considerations, we concluded the stochas- 
tic forces make the orbital parameters undergo a random walk 
that is dependent on the force model primarily through the dif- 
fusion coefficient D ^ 2(^Ff) Tc- 

For planets under the gravitational influence of a protoplan- 
etary disk, the natural scale for the force per unit mass compo- 
nents, Fi, is i^o(r) = TiQYi{ r)/2, where S is the ch aracteristic 
disk surface density (see eg. Papaloizou & Terquem 2006). We 
comment that Fq (r) is the gravitational force per unit mass due 
to a small circular disk patch of radius at a distance \/2 
from its centre assuming that all its mass is concentrated there. 
The result is independent of . The natural correlation time Tc 
is the inverse of the orbital angular frequency Tc_o = il^^. To set 
the natural scale for D, we adopt a minimum mass solar nebula 
model (MMSN, see lWeidenschill ing 1977) with 



I](r) = 4200- 



cm2 Vl AU/ 



-3/2 



(63) 



This provides a natural scale for 13 as a function of the local disk 
radius and the central stellar mass through 



Do = 2CFoVe,o = 25 C 



(l Au) 



(64) 



where C is a dimensionless constant. 

There are several very uncertain factors which contribute 
to determining an appropriate value of the dimensionless con- 
stant C: The density fluctuations f ound in MRI simulations are 
typically Sp/p = ST,/T, w 0.1 (eg. lNelsonir2005h . The presence 
of a dead zone in the mid plane regions of the disk, where the 
MRI is not active, has been found to cause reductions in the mag- 
nitude of Fq by one order of magnitude or more, as compared to 
active ca ses dOishi et al.l[2p07.) . Massive planets open a gap in 
the disk. lOishi et aL ( 20071) found that most of the contribution 
to the stochastic force comes from density fluctuations within a 
distance of one scaleheight from the planet. When a gap forms, 
this region is cleared of material leading to the expectation of 
a substantial decrease in the the magnitude of turbulent density 



fluctuations. Consequently Fq should be reduced on account of a 
lower ambient surface density. A factor of -,7; seems reasonable 
although it might be even smaller ( de Val-Borro et al. 2006). The 
correlation time Tc is actually found to be appr oximately . 5 ^ ^ 
jNelson & Papaloizoul2004l lOishi et al.ll2007l) . 

If it is appropriate to include reduction factors to account for 
all of the above effects, one finds C ~ 5 ■ 10~^ and we expect a 
natural scale for the diffusion coefficient to be specified through 



Dn 



10- 



(1 Au) 



-3/2 / 



IMf. 



-1/2 2 
' cm"^ 



(65) 



The same value of D may be equivalently scaled to the orbital 
parameters of the planets without reference to the disk by writing 



Dn = 3.5 • W 



-5/2 / M. 



1 AU 



I Mr. 



3/2 



cm 



(66) 



Thus a value Dq = 10~^ in cgs units corresponds to a ratio of 
the root mean square stochastic force component to that due to 
the central star of about ~ 10~^ for a central solar mass at 1 AU. 
It is a simple matter to scale to other locations. 

Of course we emphasize that the value of this quantity is very 
uncertain, a situation that is exacerbated by its proportionality to 
the square of the magnitude of the stochastic force per unit mass. 
For this reason we perform simulations for a range of D covering 
many orders of magnitude. 

3.2. Numerical Implementation 

The procedure we implemented, uses a discrete first order 
Markov process to generate a correlated noise that is continu- 
ous and added as an additional force. The Markov process is a 
statistical process which is defined by two parameters, the root 
mean square of the amplitude and the correlation time Tc dKasdinl 
1995). It has a zero mean value and has no memory. This has the 
advantage that previous values do not need to be stored. The 
autocorrelation function decays exponentially and thus mimics 
the autocor relation function measured in MHD simulations by 
lOishi et al.l(l2007l) . 

3.3. Stochastic Forces Acting on a Single Planet 

We first investigate the long term effect of stochastic forces on a 
single isolated planet. The initial orbital parameters were taken 
to be the the observed parameters of GJ876 b (see table [Til and 
the central star had a mass of 0.38 Mq. In this simulation, we 
use the reduced diffusion coefficients 



D = 8.2- 10" 



D = 8.2- 10" 



and 



(67) 
(68) 



which can be represented by a correlation time of half the or- 
bital period and a specific force with root mean square values 
^{F^) = 4.05-10"6cm/s2 and 



' cm/s 



%F2y ^ 4.05-10- 
respectively (see also equation (|66l l above). 

The resulting random walks undergone by e, a and ru are 
plotted in figure [T] for six different realisations for each of the 
two diffusion parameters. The spreading rates can be estimated 
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Fig, 1. Time evolution of the longitude of periastron tu, the ec- 
centricity e and the semi major axis a for a single planet. The 
initial orbital parameters of the planet were taken to be those of 
GJ876 b (see table[T]i. Six different realisations starting from the 
same initial conditions are shown in each panel. The diffusion 
coefficients are = 8.2 • lO^^cm^/s^ for the upper three panels 
and D = 8.2 • 10^'^cm^/s'^ for the lower three respectively. The 
solid lines correspond to the analytic predictions for the amount 
of spreading (see text). 

from equations ( |46] | - ( |49] l. These analytic predictions are plotted 
as solid lines. Clearly the numerical model is in broad agreement 
with the random walk description with a spreading that scales 
with D as expected. We have performed simulations for a variety 
of D and get results that are fully consistent with this in all cases. 

3.4. Illustration of the Modes of Libration in a Two Planet 
System 

We now go on to consider two planet systems. As an illustrative 
example we consider the GJ876 system (see table [T]i. Note that 
we can easily scale all physical quantities and extend the discus- 
sion to other systems (see section|6]for a discussion). We begin 
by considering the evolution of the system without stochastic 
forcing in order to characterize the modes of libration of the res- 
onant angles and other orbital parameters. In particular we iden- 
tify the fast and slow modes discussed in section |2".2.2| 

The time evolution of the resonant angles and the eccentric- 
ities is plotted in figure |2] Clearly visible are the slow and fast 
oscillation modes. The fast mode, which is seen to have a period 
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Fig. 2. Time evolution of the resonant angles and the eccentric- 
ities in the system GJ876 without turbulent forcing. The domi- 
nance of the fast mode with period '--^ 1.4 years in the oscilla- 
tions of 4>i, and the dominance of the slow mode with period 
~ 10 years in ( can be clearly seen. 

of about 1.4 years, dominates the librations of 62, and (pi while 
also being present in those of 02- On the other hand the slow 
mode, which is seen to have a period of about 10 years, domi- 
nates the librations of C while also being present in those of ei 
and (j)2- 

We emphasize the fact that the eccentricities of the two plan- 
ets participate in the librations and so are not constant. In partic- 
ular, the eccentricity of GJ876 b oscillates around a mean value 
of 0.03 with an amplitude Ae w 0.01 — 0.02. This behaviour 
involving the attainment of smaller values of the eccentricity 
has important consequences for stochastic evolution as discussed 
above (see section |2". 3. 3l l and see also below. 

We remark that similar behaviour occurs for all the systems 
we have studied, these have a wide range of eccentricities and 
planet masses. Indeed we note that the period scalings of the 
fast and slow mode periods with the planet masses are provided 
by equations dSTl i and (l33T l respectively. These indicate, as con- 
firmed by our simulations, that if both masses are reduced by a 
factor A then the period of the fast mode scales as VA and the 
period of the slow mode scales as A. 

4. Two Planets with Stochastic Forcing 

We now consider systems of two planets with stochastic forcing. 
For simplicity we begin by applying forcing only to the outer 
planet. We have found that adding the same form of forcing 
to the inner planet tends to speed up the evolution by approx- 
imately a factor of two without changing qualitative details. For 
illustrative purposes we again start with the GJ876 system and 
consider two diffusion coefficients D = 8.2 • lO^^cm^/s"^ and 
D = 8.2 • 10^'*cm'^/s'^. We remark that all of our simulations 
were run with constant values of D = 2{Fi)'^Tc. This was done 
maintaining the parameters {Fi) and Tc to be constant, with Tc 
being determined for the initial location of the outer planet. For 
the cases considered here, there is little orbital migration so this 
is not a significant feature. Note also that different realizations of 
the system are likely to become unstable and scatter for higher 
values of D (see below). 

The time evolution of the eccentricities is plotted in figure 
[3] We see fast oscillations superimposed on a random walk. The 
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Fig. 3. Time evolution of the eccentricities in the GJ876 sys- 
tem with turbulent forcing included. The diffusion coefficient 
is D = 8.2 • 10~^cm^/s'^ for the two upper plots and D ~ 
8.2 • 10^'*cm^/s^ for the two lower plots. The first (uppermost) 
and the third plots each show a single run. The second and fourth 
(lowermost) plots show the time averaged eccentricity for four 
runs. The averaging interval is 1000 years. 



amplitude of the oscilations, as well as the mean value, changes 
with time. 

Our simple analytic model, assumes slowly changing back- 
ground eccentricities and semi-major axes and accordingly does 
not incorporate the oscillations of the eccentricity due to the res- 
onant interaction of the planets. In order to make a comparison, 
we perform a time average over many periods to get smoothed 
quantities whose behaviour we can compare with that expected 
from equations ( |46] l - (|49] l. When this procedure is followed, the 
evolution is in reasonable accord with that expected from the an- 
alytic model provided that allowance is made for the importance 
of small values of ei in determining the growth of the libration 
amplitude of the angle between the apsidal lines of the two plan- 
ets (see equation (|49]l). The presence of this feature results in the 
behaviour of the libration amplitude being more complex than 
that implied by a process governed by a simple diffusive random 
walk. 



4. 1. Disruption of a Resonance in Stages 

Systems with mean motion commensurabilities can be in many 
different configurations. Here we describe the important evolu- 
tionary stages as they appear in a stochastically forced system 
that starts from a situation in which all the resonant angles show 
small amplitude libration. For example observations of GJ876 
suggest that the system is currently in such a state with the ratio 
of the orbital periods Pi / P2 oscillating about a mean value of 2. 

4.1 .1 . Attainment of Circulation of tine Angle Between the 
Apsidal Lines 

When the initial eccentricity of the outer planet, ei, is small the 
excitation of the libration amplitudes of the resonant angles read- 
ily brings about a situation where ei attains very small values. 
This causes the periastron difference ( to undergo large oscil- 
lations and eventually circulation (see equation (|49]l). Should 
stochastic forces cause ei to reach zero, on account of the co- 
ordinate singularity rui becomes undefined. Subsequently very 
small perturbations are able to produce a small eccentricity with 
undergoing large amplitude librations or circulation (see be- 
low). But note that as also mentioned in section im this occurs 
without a large physical perturbation to the system. Thus the oc- 
curence of this phenomenon does not imply the system ceases to 
be in a commensurability. 

In this context we note that the eccentricity of GJ876 b ini- 
tially is such that ei ^ 0.03 with values ei ^ 0.01 often being 
attained during libration cycles. Thus only a small change may 
cause the above situation to occur. In all cases we have consid- 
ered, we find that ( enters circulation prior to the fast angle 0i , 
which may remain Ubrating until that too is driven into circula- 
tion. 

4.1 .2. Attainment of Circulation of the Fast Angle 

Both before and after ( enters circulation, stochastic forcing acts 
to increase the libration amplitude of the fast mode. This mode 
dominates both the librations of the resonant angle (pi and the 
semi major axes. Eventually ipi starts circulation. Shortly after- 
wards commensurability is lost and P1/P2 starts to undergo a 
random walk with a centre that drifts away from 2. Note that it 
is possible for some realizations to re-enter commensurability. 
For systems with the masses of the observed GJ876 system, the 
most likely outcome is a scattering event that causes complete 
disruption of the system. 

4.1.3. A Numerical Illustration 

In order to illustrate the evolutionary sequence described above 
we plot results for two realisations of the evolution of the GJ876 
system in figure |4] For these runs we adopted the diffusion co- 
efficient D — 0.42cm^/s'^. In this context we note that reducing 
D increases the evolutionary time which has been found, both 
analyticaly and numerically to be cx 1 /D (see below). 

The times at which the transition from libration to circula- 
tion occurs for both the slow and fast angles are indicated by 
the vertical lines in figure |4] Several of the features discussed 
above and in section |2.3.3| can be seen in figure |4] In particular 
the tendency for the occurence of very small values of ei to be 
associated with transitions to circulation of C, can be seen for the 
realisation plotted in the lower panels at around t = 80 years 
and t — 500 years. Such episodes always seem to occur when 
the libration amplitude of the fast angle 0i is relatively boosted. 
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Fig. 4. Time evolution of the resonant angles, the period ra- 
tio P2/P1 and the eccentricity, ei, in the GJ876 system with 
stochastic forcing corresponding io D = 0.42cm^/s''. The ver- 
tical lines indicate when the angles enter circulation for a pro- 
longed period. The realisation illustrated in the lower panel scat- 
ters shortly after 0i goes into circularization. 



indicating that this plays a role in boosting the slow angle. If the 
period of time for which ei attains small values is small and 0i 
recovers small amplitude librations, the slow angle returns to cir- 
culation. Thus the attainment of long period circulation for the 
slow angle is related to the diffusion time for 0i. We also see 
from figure |4] that the angle 02 7 which has a large contribution 
from the slow mode, behaves in the same way as C as far as li- 
bration/circularization is concerned. We have verified by consid- 
ering the results from the simulations of GJ876 LM HE, which 
started with a larger value of ei, that, as expected, the attain- 
ment of circulation of the slow angle takes relatively longer in 
this case, the time approaching more closely the time when 0i 
attains circulation. Also as expected, the time when 01 attains 
circulation is not affected by the change in ei. 

4.2. Dependence on the Diffusion Coefficient 

We now consider the stability of the systems listed in table [1] as 
a function of D. These systems have a variety of masses and 
orbital eccentricities. In particular, in view of the complex inter- 
action between the resonant angles discussed above we wish to 



investigate whether the mean amplitude growth at a given time 
is indeed proportional to Z?. As also mentioned above, the value 
of the diffusion parameter D, that should be adopted, is very 
uncertain. We have therefore considered values of D ranging 
over five orders of magnitude. However, the correlation time Tc 
is always taken to be given by Tc = 0.5il^^ while the RMS 
value of the force is changed. In order to determine the "life- 
time" of a resonant angle, we monitor whether it is librating or 
circulating. Numerically, libration is defined to cease at the first 
time the angle is seen to reach absolute values largen than 2. We 
note that the angle can in general regain small values afterwards. 
However, this is a transient effect and changes the lifetime by no 
more than a factor of ^ 2 in all our simulations. 

In this context we consider the fast angle 0i and the slow 
angle C, though, as we saw above, the latter can be replaced by 
02 which exhibits the same behaviour. As 0i is the last to start 
circulating the resonance is defined to be broken at that point. 

Equations ( l60b and ( |59] l estimate the spreading of the res- 
onant angles as a function of time. We determine the times to 
attain circulation as the times to attain (A0i)^ = 4, i = 1,2 as- 
suming the initial values for the orbital elements. We plot both 
the numerical and analytical results in figure|5] To remove statis- 
tical fluctuations and obtain a mean spreading time, the numeri- 
cal values for a particular value of D were obtained by averaging 
over 60 realizations. 

From figure|5]it is apparent that the evolutionary times scale 
is oc 1/D for D varying by many orders of magnitude. The ana- 
lytic estimates for the libration survival times of 01 , dominated 
by the fast mode, obtained using equation ( l60l l adopting the ini- 
tial orbital elements and with the fast libration frequency de- 
termined from the simulations, are plotted in the upper panel of 
figure|5] These are in good agreement with the numerical results. 
However, the estimate for 02 based on equation ( |59] l using the 
initial value of ei overestimates the lifetime by a factor of at least 
^ 40 (see solid line in the lower plot of figure |5]). Futhermore 
(|59] l has no dependence on planet mass which is in clear con- 
flict with the numerical results. This, as discussed above is due 
to the temporary attainment of small eccentricities. This causes 
the disruption of the libration of C, earlier than would be pre- 
dicted assuming ei is constant. In fact this disruption occurs at 
times that can be up to a factor of 10 times shorter than those re- 
quired to disrupt 01. We estimate the lifetime of librations of 02 
by calculating the time 0i needs to reach values close to 1 (see 
dashed lines in the lower plot of figure |5]). As explained above 
in section 12.3.31 large amplitude variations of 0i are expected 
to couple to and excite the slow mode. Variations induced in the 
eccentricity ei allow this to reach zero and we lose libration of 
02 . These simple estimates are in good agreement with the nu- 
merical results and accordingly support the idea that 0i and 02 
are indeed coupled in the non linear regime. For low mass plan- 
ets this simple analytic prediction underestimates the lifetime of 
02 by a factor of ~ 2, suggesting that the coupling between the 
two modes depends slightly on the planet masses. 

To confirm our understanding of these processes, we have 
repeated the calculations with systems that are in the same state 
as the system discussed above but have higher eccentricities (see 
GJ876 LM HE and GJ876 SE HE in table[T]i. The Ufetime of the 
librating state of 0i is unchanged, as this does not depend signif- 
icantly on the eccentricity. However, the lifetime of the librating 
state of 02 is a factor of 2 — 3 longer. This is due to the fact that 
a larger excitation of 02 is needed to make ei reach small values 
in these simulations. 

The random walk description breaks down completely when 
the anticipated disruption time becomes shorter than the libra- 
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Fig. 5. Average time until circulation of the resonant angles 
01 (top) and (/)2 (bottom) in the GJ876 (indicated with +) , 
GJ876LM (indicated with x), GJ876SE (indicated with *) and 
GJ876E (indicated with □) systems as a function of the stochas- 
tic forcing diffusion coefficient D. In the case of the upper panel, 
the analytic curves explained in the text, are from top to bot- 
tom for the GJ876, GJ876LM, GJ876SE and GJ876E systems 
respectively. In the lower panel the analytic curves, as explained 
in the text, are from top to bottom for the GJ876 (solid curve and 
top dashed curve), GJ876LM, GJ876SE and GJ876E systems 
respectively. 



tion period. This is expected to happen for small enough masses 
as can be verified from figure|5] It is because the disruption time 
decreases linearly with the planet masses while the libration pe- 
riod increases as the square root of the planet masses. Then we 
cannot average over many libration periods. This situation is ap- 
parent in figure |5] for very short disruption times of the order 
100-1000 years where survival times cease to vary linearly with 
l/D. 



5. Formation of HD1 28311 

We here discuss the application of the ideas discussed above to 
understanding the orbital configuration of the planetary system 
HD128311. This is (with 99% confidence) in a 2:1 mean mo- 
tion resonance with the angle 4>i librating an d the angle C, cir- 
culating so that there is no apsidal corotation dVogt et al. 2005!). 
However, the orbital par ameters are not w ell constrained. The 
orgininal Keplerian fit bv lVogt et al.l (l2005b is such that the plan- 



ets undergo a close encounter after only 2000 years. The values 
in table [T]have been obtained from a fit to the observational data 
that includes interactions between the planets. The values quoted 
have large error bars. For example the eccentricities ei and 62 
have an uncertainty of 33% and 21%, respectively. Although the 
best fit doesn't manifest apsidal corotation, the system could un- 
dergo large amplitude librations and still be stable. 

According to Lee & Peale (2002) the planets should have ap- 
sidal corotation if the commensurability was formed by the two 
planets undergoing sufficently slow convergent inward migra- 
tion, and if they then both migrated inwards sig nificantly while 
maintaining the commens urability (see also Snellgrove et alj 
2001). Accordinglv Sandor & KlevI (l2006l) suggested a possible 
formation scenario with inward migration as described above, 
but with an additional perturbing event, such as a close encounter 
with an additional planet occuring after the halting of the inward 
migration. This perturbation is needed to alter the behaviour of 
(, so that it undergoes circulation rather than libration and thus 
producing orbital parameters similar to the observed ones. 

We showed above that stochastic forcing possibly resulting 
from turbulence driven by the MRl readily produces systems 
with commensurabilities without apsidal corotation if the eccen- 
tricities are not too large. This suggests that a scenario that froms 
the commensurability through disk induced inward convergent 
migration might readily produce commensurable systems with- 
out apsidal corotation if stochastic forcing is included. Such sce- 
narios are investigated in this section. 

We present simulations of the formation of HD 1283 11 that 
include migration and stochastic forcing due to turbulence but 
do not invoke special perturbation events involving additional 
planets. We find that model systems with orbital parameters re- 
sembling the observed ones are readi ly produced that giv e bet- 
ter matches than so far provided by S andor & Kiev! ( 12006 ). The 
planets in this system are of the order of several Jupiter masses 
and the eccentricity of one planet can get very small during one 
libration period. The observed orbital parameters are given in ta- 
ble[T]and plotted on the right hand side of figure s|6]|7] and [8] The 
mass of the star is 0.8 M0. 

5. 1. Migration Forces 

We incorporate the non-conservative forces excerted by a pro- 
tostellar d isk that l ead to inward migration by following the ap- 
proach of iLee & Peale (2001) and Snellgrove et al. (2001). In 
this the migration process is characterized by migration and ec- 
centricity damping timescales, Ta — a /a and = e/e. We 
also define the ratio of the above timescales to be if = \Ta \ /Tg. 
This ratio determines the eccentricities in the state of self sim- 
ilar inward migration of the commensurable system that is at- 
tained after large times. We verify this result for two different 
values of K (see below) . The procedure is now widely used (e.g. 
Sandor & Kley 2006; Crida et al.1 [20081) and we have checked 
our code against the results of their work. 

In this work we allow the planets to form a commensurabil- 
ity through convergent inward migration but stochastic forcing 
is included, the disk is then removed through having its surface 
density reduced to zero on a 2000 year timescale. This simul- 
taneously reduces both the migration forces and the stochastic 
forcing to zero. This procedure is very similar to that adopted by 
the above authors but we have included stochastic forcing and 
removed the disk on longer time scales. 

The stochastic forces have to have the right balance with re- 
spect to the migration rate. We have found that inward migra- 
tion imparts stability to the resonant system. If the migration rate 
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Fig. 6. The plots on the left hand side show a possible forma- 
tion scenario of HD128311. We plot the observed system on the 
right hand side as a comparison (see table[T]and text). The plots 
show the resonant angles 01, (/)2, the eccentricities ei, 62 and 
period ratio Pi / P2 for formation scenarios including turbulence 
and migration. Resonace capturing occurs after 2000 years. 
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Fig. 7. The plots on the left hand side show another possible 
formation scenario of HD128311. Again, we plot the observed 
system on the right hand side as a comparison (see table [T] and 
text. The migration timescales in this run are Ta.i — 16000, 
Tq 2 — —40000 and K = 5.5. Note that the eccentricities are 
larger compared to figure |6]because K is smaller 



is too fast relative to the stochastic forcing the migration keeps 
down the libration amplitudes and we do not get circulation. On 
the other hand large eccentricity damping favours broken apsi- 
dal corotation. This might look counterintuitive at first sight, but 
remember that the diffusion of ( depends on 2- Due to the 
stochastic nature of the problem, it is hard to present a continuum 
of solutions so we restrict ourselves to the discussion of three 
representative examples. However, we comment that we are able 
to obtain similar end states for a wide range of migration param- 
eters. Before presenting these examples we briefly indicate how 
the librations can be stabilized by the migration process. 

As above, the parameters (Fj) and Tc are kept constant, so 
maintaining D constant, for the duration of the simulation, with 
Tc being determined for the initial location of the outer planet. 
As discussed above, these values may scale with the radius of 
the planet and we thus expect them to change during migration. 
However, the semi major axis of the outer planet changes only 
by ~ 30 % during the simulation consequently we ignore this 
effect. 




5000 10000 15000 20000 25000 30000 35000 40000 4000 8000 



1 (years) t (years) 

Fig. 8. These plots show another possible formation scenario of 
HD12831 1. The damping parameters are the same as in figure|2] 



5.2. Adiabatic Invariance and the Stabilization of 
Librations tlirougli Damping 

Both the fast angle (pi and the slow angle 02 obey stochastically 
forced oscillator equations. Although we do not have a simple 
expressions for the oscillation frequencies we know that for fixed 
eccentricity and planet masses they both scale as the orbital fre- 
quency. Til. When the planets migrate inwards together, this in- 
creases slowly with time. Accordingly the theory of adiabatic 
invariance indicates that the Ubration amplitude should scale as 
Ui (lPealelll97^ . If there were no stochastic forcing or other 
effects, the libration amplitude would decrease on twice the mi- 
gration time scale 2ra once steady eccentricities were formed. 
This provides only a small effect given the typical amount of 
migration ^ 30% in our simulations. 

A comparison of the behaviour of a model of the GJ876 sys- 
tem undergoing convergent orbital migration resulting in the fo- 
mation of a commensurabil ty both with and with out eccentricity 
damping was undertaken bv lLee & Peald (1200 lb . Small libration 



amplitudes at a similar level are attained in both models. The 
case without eccentricity damping attains higher eccentricities 
and somewhat larger libration amplitudes as the evolution con- 
tinues whereas a low steady libration amplitude is maintained 
in the case with eccentricity damping. The above situation is 
consistent with the initial reduction in libration amplitude that 
occurs prior to the attainment of steady eccentricities as not be- 
ing due to eccentricity damping but the operation of adiabatic 
invariance CPeale 1976} . 

However, a study of simulations of the GJ876 system with- 
out stochastic forcing reveals that eccentricty damping does play 
a role in damping the slow mode (associated with apsidal coro- 
tation) on the eccentricity damping timescale, whereas the fast 
mode appears to be only marginally affected at least at small 
amplitudes. Thus the residual librations in these cases consist 
mainly if not entirely of the fast mode. Therefore the inclusion 
of eccentricity damping does impart some stability to the com- 
mensurability through control of the slow mode and thus its pos- 
sible interactions with the fast mode. As the slow mode involves 
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and ei, (see figure |2]i, it is not surprising that this is affected 
by eccentricity damping. 

Note too that because the migration process itself causes li- 
bration amplitude reduction on formation of the commensurabil- 
ity some stability is also implied towards some types of distur- 
bance that would take the system away from resonance. 

5.3. Model 1 

The planets are initially in circular orbits at radii ai = 4 AU and 
a2 = 2 AU. Stochastic forcing is applied to the outer planet only 
with the diffusion coefficient D — 6.4- 10^'^cm^/s'^. Note that 
although this value is 640 times larger than the scaling given by 
equation ( l65T l. corresponding to a force that is 25 times larger 
than the simple estimate given in section 13.11 corresponding to 
smaller reduction factors resulting from a gap or dead zone, it 
is too small for the angle t/fi to be brought to circulation during 
our runs. The results given above and summarized in figure |5] 
and other tests, indicate that similar results would be obtained if 
D is reduced, but on a longer time scale cx 1/D, provided the 
migration rate is also appropriately ultimately reduced so that 
the system can survive for long enough to enable C, to be driven 
into circulation. 

The outer planet is made to migrate inwards on a timescale 
— 8000 years. The inner planet migrated slowly outwards 
on a timescaly of — —20000 years. This results in con- 
vergent migration. For both planets we use an eccentricity to 
semi major axis damping ratio of K = 8. The resulting evo- 
lut ionary timescales are s ignificantly larger than those used by 
eg. ISandor&KlevI (120061) and more easily justified by hydrody- 
namical simulations. 

The time evolution is shown in the left plot of figure|6] After 
resonance capturing all angles are either initially librating or on 
the border between libration and circulation. The slow mode 
has a period of ~ 700 years, whereas the slow mode period is 
20 times shorter 

It can be seen from figure |6] and other related figures be- 
low, that while migration continues libration ampliudes tend to 
be controlled apart from when ei becomes either zero or very 
close to zero. Then, due to stochastic forcing, 02 and C start cir- 
culating. Subsequently libration is recovered over time intervals 
for which ei does not attain very small values but eventually ad- 
ditional stochastic forcing together with the repeated attainment 
of small values for ei causes (p2 to remain circulating for the re- 
mainder of the simulation. After 13000 years, both the forces due 
to migration and turbulence are reduced smoothly on a timescale 
of 2000 years. The result is a stable configuration that resembles 
the observed system very well. 

5.4. Model 2 

For this simulation, we lengthened the migration timescales by 
a factor of 2 to show that the results are generic. We also de- 
creased the value of K to 5.5. This results in larger eccentric- 
ities and the final state better resembled our representation of 
the observed system. However, it should be kept in mind that 
the eccentricities are not well constrained by the observations. 
All other parameters are the same as for model 1 . The time evo- 
lution is plotted in figure [T] The resonant angles librate imme- 
diately after the capture into resonance as predicted for suffi- 
ciently slow migration. However, in the same way as for model 
1 described above, stochastic forces make 02 circulate soon af- 
terwards. Once the migration forces and stochastic forces are 



removed between 20000 and 22000 years, the system stays in a 
stable configuration with no apsidal corotation. 

5.5. Models 

All parameters are the same as for model 2. Accordingly model 
3 represents another statistical realization of that case. The time 
evolution is plotted in figure|8] The final state is very close to the 
boundary between libration and circulation of (. As discussed 
above, large uncertainties in the orbital parameters means that 
the actual system could be in such a state. 

Due to stochastic forces, which may result from local turbu- 
lence, we are able to generate a broad spectrum of model sys- 
tems. Some of them undergo a strong scattering during the mi- 
gration process. However, the surving systems resemble the ob- 
served system very well. Although the orbital parameters are not 
well constrained yet, we showed that values in the right range are 
naturally produced by a turbulent disk. 

6. Discussion 

In this paper, we have presented a self consistent analytic model 
applicable to either a single planet or two planets in a mean 
motion resonance subject to external stochastic forcing. The 
stochastic forces could result from MRl driven turbulence within 
the protoplanetary disk but our treatment is equally applicable to 
any other source. 

We considered the evolution of a stochasticly forced two 
planet system that is initially deep inside a MMR (ie. the 
two independent resonant angles librate with small amplitude). 
Stochastic forces cause libration amplitudes to increase in the 
mean with time until all resonant angles are driven into circula- 
tion at which point the commensurability is lost. Often a strong 
scattering occurs soon afterwards for sytems composed of plan- 
ets in the Jovian mass range. 

We isolated a fast libration mode which is associated with 
oscillations of the semi-major axes and a slow libration mode 
which is mostly associated with the motion of the angle between 
the apsidal lines of the two planets. These modes respond dif- 
ferently to stochastic forcing, the slow mode being more readily 
converted to circulation than the fast mode. This slow mode is 
sensitive to the attainment of small eccentricities which cause 
rapid transitions between libration and circulation. The ampli- 
tude of the fast mode grows more regularly in the mean, with 
the square of the libration amplitude in most cases increasing 
linearly with time and being proportional to the diffusion coeffi- 
cient D. Of course this discussion is simplified and there are lim- 
itations. For example if the total mass of the system is reduced, 
the disruption time eventually becomes comparable to the libra- 
tion period. In that case the averaging process that we used in the 
derivation is no longer valid and the lifetime no longer scales as 
1/D. 

The analytic model was compared with numerical simula- 
tions which incorporated stochastic forces. Those forces, param- 
eterized by the mean square values of each component in cylin- 
drical ccoordinates and the autocorrelation time, were applied in 
a continuous manner giving results that could be directly com- 
pared with the analytic model. The simulations were in broad 
agreement with analytic predictions and we presented illustra- 
tive examples of the disruption process. We performed the sim- 
ulations for a large range of diffusion parameters, planet masses 
and initial eccentricities to verify the scaling law for the com- 
mensurability disruption time summarized hereafter. 
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To summarize our results, recalling that the slow angle is 
driven into circulation before the fast angle, so that the ultimate 
lifetime is determined by the time taken for the fast angle to 
achieve circulation, we determine the lifetime, tf, using equa- 
tion (|60]l setting t ^ tf, (AQ)^ {^Q)o 4, together with 
7/ = p = 1, so obtaining 



/ 



36D 



(69) 



We showed above that this gives good agreement with our nu- 
merical results. Using D = 2{F^)tc, we can express this result 
in terms of the relative magnitude of the stochastic forcing in the 
form 



tf = 2.4 X 10 



8.5uJify/qGJ 



V 8ni 



qoj 



-Pi, 



(70) 



where Pi is the orbital period of the outer planet. Here the first 
quantity in brackets represents the ratio of the square of the cen- 
tral force per unit mass to the mean square stochastic force per 
unit mass acting on the outer planet. The other quantities in 
brackets, scaled to the GJ876 system are expected to be unity, 
while the last factor q/qa.j is the ratio of the total mass ratio of 
the system to the same quantity for GJ876. Here it is assumed 
that the two planets in the system have comparable masses. 

From ( |70l ) we see that a non migrating system such as GJ876 
could survive in resonance for tf 10^ years if the stochas- 
tic force amplitude is ^ 10~^ times the central force. This ex- 
pression enables scaling to other systems at other disk loca- 
tions for other stochastic forcing amplitudes. Inference of sur- 
vival probabilities for particular systems depends on many un- 
certain aspects, such as the protoplanetary disk model and the 
strength of MRI turbulence. However, the mass ratio dependence 
in equation jTOl ) indicates that survival is favored for more mas- 
sive systems. At the present time the number of observed reso- 
nant systems is too small for definitive conclusions to be made. 
However, the fact that several systems exhibiting commensu- 
rabilities have been observed indicates that resonances are not 
always completely disrupted by stochastic forces due to turbu- 
lence, but rather may be modified as in our study of HD128311. 

The HD128311 system is such that the fast mode librates 
with the slow mode being near the borderline between libra- 
tion and circulation. We found that such a configuration was 
readily produced in a scenario in which the commensurability 
was formed through a temporary period of convergent migration 
with the addition of stochastic forcing. During a migration phase 
moderate adiabatic invariance applied to the libration modes to- 
gether with eccentricity damping leads to increased stabilization 
and a longer lifetime for the resonance. However, as discussed 
above, the time evolution of the eccentricity and in particular 
the attainment of small values plays an important role in causing 
the slow mode to circulate corresponding to the loss of apsidal 
corotation. Thus we expect that a large eccentricity damping rate 
does not necessarily stabilize the apsidal corotation of the sys- 
tem. Additional simulations have shown that this is lost more 
easily for large damping rates (K 3> 10). On the other hand it 
should be noted that this has less significance when one of the 
orbits becomes nearly circular 

Further observations of extrasolar planetary systems leading 
to better statistics may lead to an improved situation for assess- 
ing the role of stochastic forcing in future. 
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